## 

rm(list = ls())

library(tidyverse)
library(pbapply)
library(rdrobust)
library(fixest)
library(broom)

## Def covars

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')

## Def helper function

make_binary <- function (v, method = "median", q = 0.75) 
{
  if (!(method %in% c("mean", "median", "quantile"))) {
    (break)("Unknown method")
  }
  if (!method == "quantile") {
    cutoff = ifelse(method == "mean", mean(v, na.rm = T), 
                    median(v, na.rm = T))
    b <- ifelse(v > cutoff, 1, 0)
  }
  else {
    cutoff = quantile(v, probs = q, na.rm = T)
    b <- ifelse(v > cutoff, 1, 0)
  }
  b
}

## Get federal election data

bt <- readRDS("data/data_federal.rds") %>%
  filter(!is.na(treated)) 

## Load income data ##

inc <- read_rds("data/county_wage_gdp_unem.rds") %>% 
  dplyr::select(ags, year, matches("wage|gdp")) %>% 
  filter(year == 2013) %>% 
  dplyr::select(-year) %>% 
  mutate(across(-ags, make_binary)) %>% 
  dplyr::rename(county_id = ags)

## Merge

bt <- bt %>% 
  mutate(county_id = substr(ags, 1 ,5)) %>% 
  left_join(inc)

## List of outcomes

outcomes <- c("turnout_party", 
              "agg_left_party", 
              "agg_center_party",
              "agg_right_party")

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## Drop the 09 obs

bt <- bt %>% 
  filter(year > 2012)

## Prep for RDD; first differences

diff_df <- pblapply(outcomes, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  mutate(county_id = substr(ags, 1, 5)) %>% 
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id,
                                 gdp_nominal_capita,
                                 wage_nominal_gross) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>% 
  mutate(east = as.numeric(substr(ags, 1, 2)) > 11)

## Make separate DFs

list_dfs <- list(diff_df %>% filter(gdp_nominal_capita == 1), 
                 diff_df %>% filter(gdp_nominal_capita == 0), 
                 diff_df %>% filter(wage_nominal_gross == 1), 
                 diff_df %>% filter(wage_nominal_gross == 0))
list_labels <- c('Above-median GDP/capita', 'Below-median GDP/capita',
                 'Above-median wages', 'Below-median wages')

## Estimate across subsets and outcomes

out_rd_opt <- pblapply(outcomes, function(o) {
  lapply(1:length(list_dfs), function(i) {
    
    df_use <- list_dfs[[i]]
    df_label <- list_labels[i]
    
    ## Estimate
    
    out <- rdrobust(y = df_use %>% pull(!!o), 
                    x = df_use$runvar, c = 0,
                    covs = df_use[, covs])
    
    ## Return this : Robust B-C SE
    out2 <- out %>% 
      tidy_rd(se_nr = 3) %>% 
      mutate(bw_bias = out$bws[2,1],
             n = sum(out$N_h),
             sample = df_label,
             outcome = o,
             sample = df_label)
    
    ## Both
    
    out2
  }) %>% reduce(rbind)
}) %>% reduce(rbind)

## Some additional transformations

out_rd_opt <- out_rd_opt %>%
  filter(outcome %in% c("turnout_party",
                        "agg_left_party")) %>% 
  mutate(mod = ifelse(str_detect(sample, 'GDP'),
                      "Moderator: GDP/capita",
                      "Moderator: Wages")) %>% 
  mutate(outcome = ifelse(str_detect(outcome, "left"),
                          "Outcome: left−wing parties",
                          "Outcome: turnout"))

# Figure A.23 : effect conditional on GDP/wages ----

pd <- position_dodge(0.4)

p1 <- ggplot(out_rd_opt, 
             aes(sample, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.65) +
  geom_point(shape = 21, fill = 'white', size = 2,
             position = pd) +
  xlab('') + ylab('RD estimate (percentage points)') + theme_bw() +
  facet_wrap(~outcome) +
  theme(legend.position = 'bottom') +
  coord_flip() +
  xlab("Sample")
p1


